{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pore-Scale Models\n",
    "\n",
    "Pore scale models are one of the more important facets of OpenPNM, but they can be a bit confusing at first, since they work 'behind-the-scenes'.  \n",
    "They offer 3 main advantages:\n",
    "\n",
    "1. A large library of pre-written models is included and they can be mixed together and their parameters edited to get a desired overall result.\n",
    "2. They allow automatic regeneration of all dependent properties when something 'upstream' is changed.\n",
    "3. The pore-scale model machinery was designed to allow easy use of custom written code for cases where a prewritten model is not available. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best way to explain their importance is via illustration.  \n",
    "\n",
    "Consider a diffusion simulation, where the diffusive conductance is defined as:\n",
    "\n",
    "$$ g_D = D_{AB}\\frac{A}{L} $$\n",
    "\n",
    "The diffusion coefficient can be predicted by the Fuller correlation:\n",
    "\n",
    "$$ D_{AB} = \\frac{10^{-3}T^{1.75}(M_1^{-1} + M_2^{-1})^{0.5}}{P[(\\Sigma_i V_{i,1})^{0.33} + (\\Sigma_i V_{i,2})^{0.33}]^2} $$\n",
    "\n",
    "Now say you want to re-run the diffusion simulation at different temperature.  This would require recalculating $D_{AB}$, followed by updating the diffusive conductance.\n",
    "\n",
    "Using pore-scale models in OpenPNM allows for simple and reliable updating of these properties, for instance within a for-loop where temperature is being varied.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Existing Pore-Scale Models\n",
    "\n",
    "The first advantage listed above is that OpenPNM includes a library of pre-written model.  In this example below we can will apply the Fuller model, without having to worry about mis-typing the equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:31.334208Z",
     "iopub.status.busy": "2021-06-24T11:25:31.332749Z",
     "iopub.status.idle": "2021-06-24T11:25:31.925227Z",
     "shell.execute_reply": "2021-06-24T11:25:31.923818Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "import openpnm as op\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "pn = op.network.Cubic(shape=[5, 5, 1], spacing=1e-4)\n",
    "geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to define the gas phase diffusivity.  We can fetch the ``fuller`` model from the ``models`` library to do this, and attach it to an empty phase object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:31.933560Z",
     "iopub.status.busy": "2021-06-24T11:25:31.932119Z",
     "iopub.status.idle": "2021-06-24T11:25:31.935015Z",
     "shell.execute_reply": "2021-06-24T11:25:31.936243Z"
    }
   },
   "outputs": [],
   "source": [
    "air = op.phases.GenericPhase(network=pn)\n",
    "f = op.models.phases.diffusivity.fuller\n",
    "air.add_model(propname='pore.diffusivity',\n",
    "              model=f,\n",
    "              MA=0.032, MB=0.028, vA=16.6, vB=17.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we had to supply the molecular weights (``MA`` and ``MB``) as well as the diffusion volumes (``vA`` and ``vB``).  This model also requires knowing the temperature and pressure, but by default it will look in 'pore.temperature' and 'pore.pressure'.\n",
    "\n",
    "Next we need to define a physics object with the diffusive conductance, which is also available in the model libary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:31.954627Z",
     "iopub.status.busy": "2021-06-24T11:25:31.948791Z",
     "iopub.status.idle": "2021-06-24T11:25:31.960890Z",
     "shell.execute_reply": "2021-06-24T11:25:31.959749Z"
    }
   },
   "outputs": [],
   "source": [
    "phys = op.physics.GenericPhysics(network=pn, phase=air, geometry=geo)\n",
    "f = op.models.physics.diffusive_conductance.ordinary_diffusion\n",
    "phys.add_model(propname='throat.diffusive_conductance',\n",
    "               model=f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly we can run the Fickian diffusion simulation to get the diffusion rate across the domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:31.972880Z",
     "iopub.status.busy": "2021-06-24T11:25:31.971479Z",
     "iopub.status.idle": "2021-06-24T11:25:32.111655Z",
     "shell.execute_reply": "2021-06-24T11:25:32.112431Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.38747949e-10]\n"
     ]
    }
   ],
   "source": [
    "fd = op.algorithms.FickianDiffusion(network=pn, phase=air)\n",
    "fd.set_value_BC(pores=pn.pores('left'), values=1)\n",
    "fd.set_value_BC(pores=pn.pores('right'), values=0)\n",
    "fd.run()\n",
    "print(fd.rate(pores=pn.pores('left')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Updating parameter on an existing model\n",
    "\n",
    "It's also easy to change parameters of a model since they are all stored on the object (``air`` in this case), meaning you don't have to reassign a new model get new parameters (although that would work).  The models and their parameters are stored under the ``models`` attribute of each object.  This is a dictionary with each model stored under the key match the ``propname`` to which is was assigned. For instance, to adjust the diffusion volumes of the Fuller model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.119402Z",
     "iopub.status.busy": "2021-06-24T11:25:32.118708Z",
     "iopub.status.idle": "2021-06-24T11:25:32.122157Z",
     "shell.execute_reply": "2021-06-24T11:25:32.121557Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diffusivity before changing parameter: 2.067547836430591e-05\n",
      "Diffusivity after: 2.0969679602405972e-05\n"
     ]
    }
   ],
   "source": [
    "print('Diffusivity before changing parameter:', air['pore.diffusivity'][0])\n",
    "air.models['pore.diffusivity']['vA'] = 15.9\n",
    "air.regenerate_models()\n",
    "print('Diffusivity after:', air['pore.diffusivity'][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Replacing an existing model with another\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say for some reason that the Fuller model is not suitable.  It's easy to go 'shopping' in the models library to retrieve a new model and replace the existing one.  In the cell below we grab the Chapman-Enskog model and simply assign it to the same ``propname`` that the Fuller model was previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.128203Z",
     "iopub.status.busy": "2021-06-24T11:25:32.127508Z",
     "iopub.status.idle": "2021-06-24T11:25:32.131548Z",
     "shell.execute_reply": "2021-06-24T11:25:32.130955Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diffusivity after: 1.5458486152825926e-05\n"
     ]
    }
   ],
   "source": [
    "f = op.models.phases.diffusivity.chapman_enskog\n",
    "air.add_model(propname='pore.diffusivity',\n",
    "              model=f, MA=0.0032, MB=0.0028, sigma_AB=3.467, omega_AB=4.1e-6)\n",
    "print('Diffusivity after:', air['pore.diffusivity'][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we don't need to explicitly call ``regenerate_models`` since this occurs automatically when a model is added.  We do however, have to regenerate ``phys`` object so it calculates the diffusive conductance with the new diffusivity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.140014Z",
     "iopub.status.busy": "2021-06-24T11:25:32.139367Z",
     "iopub.status.idle": "2021-06-24T11:25:32.161585Z",
     "shell.execute_reply": "2021-06-24T11:25:32.162114Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.03738023e-10]\n"
     ]
    }
   ],
   "source": [
    "phys.regenerate_models()\n",
    "fd.reset()\n",
    "fd.run()\n",
    "print(fd.rate(pores=pn.pores('left')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing dependent properties\n",
    "\n",
    "Now consider that you want to find the diffusion rate at higher temperature.  This requires recalculating the diffusion coefficient on ``air``, then updating the diffusive conductivity on ``phys``, and finally re-running the simulation.  Using pore-scale models this can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.167731Z",
     "iopub.status.busy": "2021-06-24T11:25:32.167066Z",
     "iopub.status.idle": "2021-06-24T11:25:32.169807Z",
     "shell.execute_reply": "2021-06-24T11:25:32.170332Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diffusivity before changing temperaure: 1.5458486152825926e-05\n",
      "Diffusivity after: 1.9929877222108533e-05\n"
     ]
    }
   ],
   "source": [
    "print('Diffusivity before changing temperaure:', air['pore.diffusivity'][0])\n",
    "air['pore.temperature'] = 353.0\n",
    "air.regenerate_models()\n",
    "print('Diffusivity after:', air['pore.diffusivity'][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the diffusivity increased with temperature as expected with the Chapman-Enskog model.  We can also propagate this change to the diffusive conductance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.177192Z",
     "iopub.status.busy": "2021-06-24T11:25:32.176560Z",
     "iopub.status.idle": "2021-06-24T11:25:32.180639Z",
     "shell.execute_reply": "2021-06-24T11:25:32.180071Z"
    }
   },
   "outputs": [],
   "source": [
    "phys.regenerate_models()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And lastly we can recalculate the diffusion rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.186817Z",
     "iopub.status.busy": "2021-06-24T11:25:32.186185Z",
     "iopub.status.idle": "2021-06-24T11:25:32.197590Z",
     "shell.execute_reply": "2021-06-24T11:25:32.198115Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.33744407e-10]\n"
     ]
    }
   ],
   "source": [
    "fd.reset()\n",
    "fd.run()\n",
    "print(fd.rate(pores=pn.pores('left')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating Custom Models\n",
    "\n",
    "Lastly, let's illustrate the ease with which a custom pore-scale model can be defined and used.  Let's create a very basic (and incorrect) model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.202961Z",
     "iopub.status.busy": "2021-06-24T11:25:32.202310Z",
     "iopub.status.idle": "2021-06-24T11:25:32.204196Z",
     "shell.execute_reply": "2021-06-24T11:25:32.204717Z"
    }
   },
   "outputs": [],
   "source": [
    "def new_diffusivity(target, A, B, \n",
    "                    temperature='pore.temperature', \n",
    "                    pressure='pore.pressure'):\n",
    "    T = target[temperature]\n",
    "    P = target[pressure]\n",
    "    DAB = A*T**3/(P*B)\n",
    "    return DAB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few key points to note in the above code.  \n",
    "\n",
    "1. Every model must accept a ``target`` argument since the ``regenerate_models`` mechanism assumes it is present.  The ``target`` is the object to which the model will be attached.  It allows for the looking up of necessary properties that should already be defined, like temperature and pressure.  Even if you don't use ``target`` within the function it is still required by the pore-scale model mechanism.  If it's presence annoys you, you can put a ``**kwargs`` at the end of the argument list to accept all arguments that you don't explicitly need.  \n",
    "2. The input parameters should not be arrays (like an Np-long list of temperature values).  Instead you should pass the dictionary key of the values on the ``target``.  This allows the model to lookup the latest values for each property when ``regenerate_models`` is called.  This also enables openpnm to store the model parameters as short strings rather than large arrays.\n",
    "3. The function should return either a scalar value or an array of Np or Nt length.  In the above case it returns a DAB value for each pore, depending on its local temperature and pressure in the pore.  However, if the ``temperature`` were set to ``'throat.temperature'`` and ``pressure`` to ``'throat.pressure'``, then the above function would return a DAB value for each throat and it could be used to calculate ``'throat.diffusivity'``.\n",
    "4. This function can be placed at the top of the script in which it is used, or it can be placed in a separate file and imported into the script with ``from my_models import new_diffusivity``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add this model to our ``air`` phase and inspect the new values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:25:32.210180Z",
     "iopub.status.busy": "2021-06-24T11:25:32.209519Z",
     "iopub.status.idle": "2021-06-24T11:25:32.211859Z",
     "shell.execute_reply": "2021-06-24T11:25:32.212390Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.06722719e-05 2.06722719e-05 2.06722719e-05 2.06722719e-05\n",
      " 2.06722719e-05 2.06722719e-05 2.06722719e-05 2.06722719e-05\n",
      " 2.06722719e-05 2.06722719e-05 2.06722719e-05 2.06722719e-05\n",
      " 2.06722719e-05 2.06722719e-05 2.06722719e-05 2.06722719e-05\n",
      " 2.06722719e-05 2.06722719e-05 2.06722719e-05 2.06722719e-05\n",
      " 2.06722719e-05 2.06722719e-05 2.06722719e-05 2.06722719e-05\n",
      " 2.06722719e-05]\n"
     ]
    }
   ],
   "source": [
    "air.add_model(propname='pore.diffusivity',\n",
    "              model=new_diffusivity,\n",
    "              A=1e-6, B=21)\n",
    "print(air['pore.diffusivity'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
